ROMA: a map-making algorithm for polarised CMB data sets 
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Abstract We present ROMA, a parallel code to produce joint optimal temperature and polarisation maps out of multidetector 
CMB observations. ROMA is a fast, accurate and robust implementation of the iterative generalised least squares approach to 
map-making. We benchmark ROMA on realistic simulated data from the last, polarisation sensitive, flight of BOOMERanG. 
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1. Introduction 

A direct consequence of the presence of anisotropics in the 
Cosmic Microwave Background (CMB) radiation is a certain 
degree of polarisation dRees 19681 [Kaiser, 1983| l. Constraining 
CMB polarisation provides valuable cosmological information, 
complementary to that encoded in temperature anisotropy and 
it significantly tightens the constraints on cosmological param- 
eters (Seljak & Zaldarriaga 1996'| ». 

Unfortunately, the polarised component of the CMB is ex- 
pected to be small compared to total intensity, making its 
measurement an experimental challenge. For this reason, un- 
til very recently the experimental eff'ort has not focused on 
the polarised component. However, the situation is quickly 
changing. A first detection of CMB polarisation, in agreement 
with theoretical predictions, has been announced by the inter- 
ferometric experiment DASI (IKovac et al 20 02 1. The WMAP 
satellite (IBennett et al . 200?) has measured the predicted cor- 
relation spectrum between CMB temperature anisotropy and 
polarisation (IBennett et al. 20031 . Many other CMB polari- 
sation experiments (DASI, CBI, CAPMAP, BICEP, QUEST, 
BOOMERanG 2K, MAXIPOL, SPORT, BarSPORT, Planck, 
etc.^) have either already taken data or are expected to do so in 
the very near future. 

It is well known that in CMB science the data reduc- 
tion process requires almost as much care as data gath- 
ering. The analysis of temperature data has now reached 
full maturity, and it has been proven successful on several 
datasets. The same is not true for polarisation. Although 
the overall analysis scheme can be borrowed from the tem- 
perature case, and many of the algorithms involved admit 
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a more or less straightforward generalisation, polarised data 
sets present peculiar aspects that call for specific treatments. 
One basic problem is that the tensor-like nature of polarisa- 
tion poses some constraints on the measurement process (see 
ICabella & Kamionkowski 20041 or |Lin & Wandelt, 2QQ4| for a 
recent review). As a consequence, not all possible instrumental 
configurations are equally advantageous from the point of view 
of data reduction ^Couchot et al. 19991 . The so called "opti- 
mal" configurations comprise several diff'erent polarisation sen- 
sitive detectors (polarimeters), placed at a convenient recipro- 
cal orientation. Usually data streams from diff'erent polarime- 
ters are jointly analysed to seek for a faint signal. The difficul- 
ties involved in this process, and their relation to the design of 
robust and efficient data mining algorithms, are seldom consid- 
ered in the literature (see IRevenu et al. 2000l for a noticeable 
exception). 

Two major steps in the data analysis process are (1) the pro- 
duction of sky maps from a set of Time Ordered Data (TOD) 
and (2) the extraction of the angular power spectrum from such 
maps. Here we focus on the first step, leaving the discussion 
about power spectrum estimation to a forthcoming paper. In 
a previous paper (Natoli et al.lOOl) we described the imple- 
mentation of an optimal map-making algorithm for the case of 
temperature only TOD collected by a one horned CMB experi- 
ment, and we showed the feasibility of the method by applying 
it to Planck and BOOMERanG simulated data. Here, we dis- 
cuss the generalisation of this algorithm to the case of TOD 
produced by an arbitrary number of polarisation sensitive de- 
tectors. The software implementation of this method, which we 
call ROMA (Roma Optimal Map-making Algorithm), allows 
fast and efficient production of optimal multidetector maps of 
CMB total intensity and polarisation. The code is currently 
being used to analyse the polarised data set from the second 
Antarctic flight of the BOOMERanG experiment (hereafter 
B2K, Mon troy et al. 2003l », which took place in January 2003. 
ROMA is part of a full, end-to-end data analysis pipeline for 



B2K, completely developed m Rome, whose details will be de- 
scribed elsewhere. 

The plan of this paper is as follows: In Sect. 2 we derive 
the least squares equations for multidetector map-making; in 
Sect 3 we describe the details of the implementation, while in 
Sect. 4 we report the application of ROMA to highly realis- 
tic B2K simulated data. Finally, in Sect. 5, we draw our main 
conclusions. 



2. Polarised Data Map-Making 

In this section we derive the Generalised Least Squares (GLS) 
equations for the polarised map-making problem, given an ar- 
bitrary number of polarimeters. With "polarimeter" we mean, 
here and in the following, a generic detector measuring total 
intensity plus a linear polarisation component. Other types of 
detectors do exist and rely on different strategies to measure 
polarisation. We prefer to focus on the linear polarimeter case 
because of its widespread adoption. In any case it would be 
straightforward to generalis of scheme once the data model is 
properly modified. 

The sky signal seen by a polarimeter can be expressed 
as dChandrasekhar 196Ql» : 



2cos2^-h t/sin2^). 



(1) 



where /, Q and U are the Stokes parameters for total intensity 
and linear polarisation, and the angle ^ identifies the polarime- 
ter orientation with respect to the chosen celestial frame. Note 
that we do not include the contribution of circular polarisation, 
associated to the Stokes parameter V. As a consequence, a V 
signal would be seen by our polarimeter as a contribution to /. 
This fact however is not a problem in CMB science because 
circular polarisation is not generated by Thomson scattering of 
unpolarised radiation. 

All three relevant Stokes parameters can be extracted by 
either combining the output of many detectors with diff'erent 
mutual orientation, or by allowing enough focal plane rotation. 
The data stream output of a given detector is a combination of 
sky signal and instrumental correlated noise: 



= \Atp (ip + Qp cos 2(pt + Up sin 20^) -h Ut 



(2) 



Here the index t labels time while p runs over the pixelized 
images of /, Q and U. The pointing matrix Atp couples the 
time and pixel domain. The information about beam smear- 
ing can, in general, be included in this matrix. We prefer, 
however, to assume that A is pointwise (/. e. having a single 
nonzero entry per row, occurring when a pixel falls into the 
line of sight) and hide the beam in the /, 2, and U map (/. e. 
solve for a pixelized, beam smeared image of the sky). This is 
only meaningful if the beam is, to good approximation, sym- 
metric with respect to boresight, and common to /, Q and U 
(see jArmitage & Wandelt 2004 for an algorithm that takes into 
account the effect of an asymmetric beam). In Eq. O above rit 
is a vector of correlated noise. 



By inserting the trigonometric runctions withm the pointing 
matrix, and considering a set of k polarimeters, one can recast 
Eq. O into a more compact formalism: 



Dt = AtpSp -h n 

where the datastreams of each polarimeter are combined: 
and the generalised pointing matrix becomes: 



(3) 
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^A\p cos2(ptA^p sm2(ptA^p , 
Similarly, the sky signal can be expressed as: 



ftp' 

Qp 



while the noise stream is: 



n. = 



Eq. (ISj defines a generic linear algebra problem, whose un- 
known parameters (the map pixel values) can be constrained 
by means of the standard GLS solution (e. g. jLupton, 1993| »: 

-1 



Sp = (a'n-^a)" a'n-^d. 



with: 



(4) 



(5) 



and (•) denotes the expectation value. The N matrix becomes 
block diagonal when assuming that the noise in different po- 
larimeters is uncorrected: 



{n\n^^) ={ni4) = (i ^ j) . 



(6) 



We defer the actual implementation of Eq. ^ to the next sec- 
tion. 

The treatment above is clearly idealistic. However, it is 
possible to incorporate into the formalism the correction for 
a nominal level of cross polarisation, one of the most com- 
mon systematic effects occurring in CMB polarimetry. Cross- 
polarisation is due to cross-talk between the two orthogonal po- 
larisation modes. That is, a polarimeter may be sensitive, with 
efficiency ;^poi, to radiation linearly polarised 90° off its natu- 
ral polarisation mode. If we assume that the cross-polarisation 
contamination is measurable (by on-ground and/or in-flight 
tests) and it is constant across the mission lifetime, the formal- 
ism expressed above can be generalised to take the effect into 



account. Ir we introduce a calibration constant C and a cross- 
polarisation factor ;^poi, the data model Eq. Q can be gener- 
alised as 

l)j = I [{ip + Qp cos 2(f>i + Up sin 201) 

+ Xpoi {ip - Qp COS 2(p\-Up sin 20^)] -h n\ 
that is 



tp 



1 -;\rpoi 
1 +;\rpoi 



(Qp cos20j-F t/^sin20j) 



+ (7) 



where we have embedded the cross-polarisation factor in the 
new calibration constant = c(l +;tpoi). The new data 
model expressed in Eq. is then solved by Eq. ^ pro- 
vided the generalised pointing matrix A is rewritten accord- 
ingly. Note also that if Xpo\ is close to unity , the map-making 
problem expressed by Eq. ^ becomes ill-conditioned since 
the matrix A^N"^A is singular if ;^poi = 1. For real-life exper- 
iments we expect a cross-polarisation level of < 10% {e. g. 
[Montroy et al. 2Q03| for B2K), not enough to hinder the search 
for a solution of Eq. (|4|). It is important to realize that the 
scheme outlined above can correct for a nominal, known cross 
polarisation level. Any uncertainty on this value, as well as on 
the overall calibration, is a potential source of systematics. 

3. Implementation 

The straightforward implementation of Eq. (|4|) for a modern 
CMB experiment would require storing and inverting a huge 
matrix (the map noise correlation matrix), a task often beyond 
the reach of present day supercomputers. On the contrary, it- 
erative methods only require matrix to vector products and ap- 
pear well suited to tackle the problem. One such scheme, pro- 
posed by Natoli et al. (2001), has been shown to be particu- 
larly convenient in the case of temperature map-making even 
for a high resolution full sky survey such as Planck. The idea 
( |Wright 1996") is to decompose the product A^N"^A ("unroll, 
convolve and rebin") to avoid computing and storing the map 
correlation matrix, and make use of a preconditioned conjugate 
gradient (PCG) solver. The key assumptions are: (1) assume 
that the beam is axisymmetric, so to keep the structure of the 
pointing matrix A as simple as possible and (2) assume that the 
noise is (piecewise, at least) stationary and that its correlation 
function decays after a time lag much shorter than the length of 
the timeline piece being processed. Under the last assumption, 
the N matrix is approximately circulant, and as such it is diag- 
onalized by a Fourier transform. This approach achieves linear 
scaling (with the number of time samples) per PCG iteration; 
if the system is preconditioned by the pixel hits counter (the 
number of times each pixel is seen), an accurate solution can 
be obtained in a few tens of iterations fsee lNatoli et for 
a discussion of algorithmic details; other implementations of 
this approach include .Pore et al. 2001. and Cantalupo et al.\ . 

3.1. The IQU approach 

Our approach is to extend the above mentioned method to po- 
larisation. Since the timestream includes contributions from 



both total intensity and linear polarisation (see bq. (|2|l) it is 
desirable to solve jointly for the temperature and polarisation 
maps, in the spirit of Eq. This approach, that we call IQU, 
is computationally more intensive than solving for / and (2, U) 
separately, but preferable for at least two reasons: first, it min- 
imises the Q and U contaminations on the / map; second it is 
more general, because it does not require any special constraint 
on the noise properties of the polarisation sensitive detectors, as 
instead would be the case when the TOD are explicitly differ- 
enced to remove the intensity component (it is easy to realize 
that in the last case the matrix N would not be circulant ex- 
cept in the special case of having identical noise properties in 
different channels). 

While the intensity (temperature) field can in principle be 
reconstructed from a single, one horned detector (up to a level 
of the same order of magnitude as the polarisation contribution 
to the timeline), the same thing is not true for polarisation, un- 
less special care is taken in ensuring that the detector observes 
all pixels of the sky under very different orientations. The lat- 
ter strategy has indeed been chosen by a few experiments (e. g. 
MAXIPOL, that uses a rotating grid, Johnson et al. 2003), but 
the majority (e. g. B2K, Planck) plan to extract the polarisa- 
tion maps by comparing channels having different mutual ori- 
entation of the plane sensitive to linear polarisation. If such a 
strategy is used, it is necessary to reduce several channels si- 
multaneously in order to obtain a polarisation map, /. e. one 
is forced to perform multi-detector map-making. Of course, 
multi-detector map-making can also be performed in the case 
of polarisation insensitive detectors, when an optimal tempera- 
ture map has to be produced out of different channels at a given 
frequency (not necessarily taken by the same experimental sys- 
tem). The IQU approach is therefore very general, as it also 
allows in a very natural way to merge data taken by different 
channels, for temperature and polarisation at the same time. 

Due to optical constraints, different detectors do not ob- 
serve exactly the same region of the sky during the mission 
lifetime, with the noticeable exception of full sky surveys. A 
polarisation map can only be extracted for pixels observed with 
sufficient variety of focal plane orientation, in order not to in- 
cur in an ill-conditioned A^N"^A matrix. Hence, it is a very 
natural choice to restrict the polarisation map to intersection 
of the sky regions surveyed by different detectors. We impose 
this constraint by flagging the time samples associated to pix- 
els outside of the common region. A more refined discrimi- 
nation would entail selecting explicitly those pixels observed 
with enough spread in the polarisation angle. For the sake of 
simplicity we choose to restrict the (2, U) maps to the intersec- 
tion region of detector coverage. The same constraint is neither 
necessary nor advantageous for the / map: it is not necessary, 
because the temperature is normally well defined outside of the 
commonly observed region (although in general with a larger 
noise contribution); it is not advantageous, because of the way 
the solver works: in fact, scans would often go continuously in 
and out of the common region; when unrolling a tentative so- 
lution map m over the timestream (/. e. when performing the 
operation Am), the latter would display time gaps if the solu- 
tion is restricted to the common region. One way to tackle this 
problem is to fill the gaps with a constrained noise realization. 



1 his approach has the disadvantage or increasing the computa- 
tional cost of the algorithm, because standard methods to pro- 
duce constrained noise realizations scale as the third power of 
the noise correlation length dHoffman & Ribak 199111 . We pre- 
fer instead, to avoid the introduction of extra gaps in the un- 
rolled timestream, to solve for / in the whole "union" region, 
and using the / tentative solution map to obtain a continuous 
timestream. 

3.2. Noise estimation 

The GLS algorithm described above requires knowledge of 
the noise correlation properties of each detector system. These 
must be estimated from the data themselves, which consist of a 
combination of noise and (often non negligible) signal. Hence, 
the problem of noise estimation. One possible solution is to 
start from a non optimal yet unbiased estimate of the signal 
(e. g. obtained by naive coadding), subtract it from the data 
and use the residual to produce an estimate of the underlying 
noise power spectrum. The latter can be used to compute a GLS 
map, employed in turn to produce a new noise estimate, and the 
scheme can then be iterated. This is basically the approach de- 
scribed in Prunet et al (2000) and Dore et al (2001), and is 
de facto identical to a scheme proposed earlier by Ferreira & 
Jaffe (2000), although this latter algorithm is cast in a Bayesian 
formalism. An alternative approach, proposed by Natoli et al. 
(2002) is to stop the noise iteration at the first step (thus sav- 
ing computational time) and increase statistical efficiency by 
performing a parametric fit on the estimated noise power spec- 
tral density. Although computationally more efficient, this lat- 
ter scheme depends somehow on the parametric assumption for 
the functional form of the noise power spectral density. For the 
purpose of this paper we do not want to loose generality, and 
therefore we opt for estimating the noise by means of the it- 
erative scheme. However, the parametric approach is a viable 
option for those cases where the noise spectral density can be 
conveniently modelled. 

In the case of multi-detector, polarisation sensitive map- 
making, it is a logical choice to use the global IQU maps as 
an estimate of the underlying signal. We will show below that 
this strategy achieves a substantial reduction of the estimated 
noise residual bias, of order the number of detector considered, 
as one would naively guess by applying a heuristic argument. 

A significant fraction of the data gathered from a real world 
experiment are badly contaminated and cannot be used for sci- 
ence extraction. Corruption of these data normally occurs be- 
cause of annoying events, either unpredictable {e. g. cosmic ray 
hits) or necessary {e. g. the ignition of a calibrator). Data loss 
can also aff'ect pointing information: in this case, the corre- 
sponding sky signal must be discarded even if usable per se. 
In any case, excerption of the aff'ected samples destroys the 
timeline continuity and, hence, noise correlation. Both things 
must be reintroduced somehow; this is normally done by in- 
serting a noise realization, constrained to reproduce the cor- 
rect statistical properties while being continuous (in a stochas- 
tic sense) at the gap boundaries (see Hoff'man & Ri bak 1991i 
[Stompor et al. 2002| l. At the same time, the map making code 
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Figure 1. The B2K focal plane (IMasi et al. 200311 . The 

145 GHz PSB's are located in the lower row. The upper de- 
tector row hosts unpolarised (spider web) bolometers coupled 
to a polarimeter. 

must be warned that such samples do not contain any useful 
signal, a requirement trivially accomplished by associating to 
all samples a timeline "flag". Flagged samples do not project 
over the sky but are assigned to a virtual pixel, whose value is 
estimated by the solver. Note that this scheme is equivalent to 
an ideally modified experiment, which observes the sky and, 
for flagged samples, a null signal (virtual) pixel. The latter is 
eventually assigned a noise variance (and covariance to real 
pixels) compatible with instrumental properties. As shown in 
the next section, we find that this approach allows for precise 
signal recovery in the noise free limit, while correctly minimis- 
ing the noise contribution in the weak signal regime. 

4. Application to a realistic dataset 

To test ROMA we simulated and reduced a full mission scan 
(about 10 days) of the 145 GHz channels of the B2K exper- 
iment, using the nominal pointing solution, which includes a 
polarisation angle for each detector horn (de Bernardis, priv. 
comm.). 

B2K ( |Montroy et al. 2003|» i s basically the follow up of 
BQQMERanG (ICril^^/^003l» featuring polarisation sensi- 
tive bolometers (PSB) and a scanning strategy optimised to 
measure the polarisation of the CMB and of Galactic fore- 
grounds. The portion of the sky observed by B2K overlaps 
almost completely with the target of the 1998 Antarctic cam- 
paign of BOOMERanG (de Bernardi s et al. 200011 . However, 
the scanning strategy of B2K diff'ers from the one implemented 
for the previous flight. The survey area of B2K is split between 
a ~ 300 square degrees region close to the Galactic plane and 
a low foreground emission area of ~ 1300 square degrees for 
CMB science. This area is in turn divided between a "shallow" 
and a "deep" region which difl'er by more than an order of mag- 
nitude in integration time per pixel. The deep region is totally 
embedded in the shallow one and covers about 120 square de- 
grees. 

In Fig. {I} we show the experimental setup of the B2K focal 
plane, with the 145 GHz PSB's in the lower row. 

In simulating the data the following scheme is employed: 
we generate a high resolution map (If 7 or NSIDE=2048 in 



HbALPix^ language) or a standard CMB sky wrth ACDM cos- 
mology. We assume a Gaussian beam of 10' FWHM for all the 
PSBs. The signal time streams are then obtained using the nom- 
inal B2K scan. The noise time streams are simulated assuming 
the following form for the noise spectral density: 

p(/) = A[i + (i/i//,n, (8) 

where fk is the knee frequency. We choose = 0.1 Hz, a = -2 
and an amplitude A corresponding to the expected white noise 
level of each of the 145 GHz PSB receivers; The variation 
of A between different receivers is within a factor of two (de 
Bernardis, priv. comm.). The minimum frequency in the TOD 
is set by the inverse of the mission life-time, while the nomi- 
nal sampling frequency of the experiment (60 Hz) determines 
the Nyquist frequency. For the sake of simplicity we assume 
that no cross-polarisation contamination is present in the simu- 
lations hereafter described. However, we have performed sev- 
eral tests by varying x^o\ up to a 20% level; they all concor- 
dantly show that the correct solution is found, with a small in- 
crease in the iteration count and a slightly higher residual noise 
level in the final maps. We fix the output map resolution to 3. '7 
(NSIDE=1024) to properly sample the 145 GHz beams. 

As a first step we perform noise estimation. We follow the 
iterative scheme explained above, and we compare the qual- 
ity of the noise spectrum recovered when using the single PSB 
(temperature only) maps as the signal baseline, as opposed to 
using the full (global) temperature map. Not surprisingly, in the 
latter case, the residual noise bias is substantially minimised, 
especially at intermediate frequencies (see Fig. O) when com- 
pared to the unbiased (by definition) estimate of the power 
spectrum obtained from the noise-only TOD used in the S-i-N 
simulation. 

In Fig. Q we display our results for a signal-only B2K 
timestream. Despite the highly realistic simulation of a com- 
plex experiment, the signal maps are recovered to high preci- 
sion, better than 1 % for most pixels in the I map, and better than 
10% for the Q and U maps. Note that the latter result is con- 
sistent, because the standard deviation of the CMB signal in 
the polarisation maps is roughly one order of magnitude lower 
than the one in the temperature map. Also note that in the deep 
region, the high number of hits per pixel allows to recover the 
polarisation pattern to much greater precision. 

The "signal plus noise" (S-hN) maps are displayed in 
Fig. (|4|). We have accurately tested that the code preserves the 
linearity of the GLS solution it implements: that is, running 
ROMA on a S-l-N timestream is equivalent to summing the out- 
put of two separate S and N runs. Strictly speaking the latter 
statement is true only if the code is allowed to iterate until the 
final solution is reached. We found in practice that the linearity 
is very well preserved after a reasonable number of iteration 
(~ 100) are completed. 

The importance of having a joint IQU solver is stressed 
by Fig. (121). Here we present, in the same fashion of the bot- 
tom row of Fig. ^ above, the diff'erence (input minus output) 
maps for a signal only case. However (see also the figure's cap- 
tion), we show here the residuals obtained when processing a 

^ http : //www . eso . org/science/healpix/ 
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Figure 2. Shown, as a function of frequency, is the percentual 
ratio between the noise power spectrum of a single detec- 
tor estimated from the 5" + N data (by means of the iterative 
procedure described in the text) and the unbiased estimate of 
the same noise power spectrum obtained from the noise-only 
timestream. The black (lower) curve is for a single detector 
case, while the blue (upper) curve is obtained with a GLS map 
from all eight detectors. Note the strong suppression of the 
residual bias in this latter case. 

single channel map, for which no polarisation solution can be 
found. As one would expect, this residual is dominated by a 
polarisation-like pattern, which the map-making code is unable 
to distinguish from the temperature signal. 

Most of the computational eff'ort required by ROMA is 
claimed by the Fourier transforms, needed to perform repeated 
convolutions with the N matrix. An efficient FFT library must 
therefore be used. Our choice falls on the publicly available 
FFTW3 library ( |Frigo & Johnson 1 998|l, which claims, in our 
case, about 80% of the total computing time. Use of the FFT 
guarantees that the code scales linearly (per PCG iteration) 
with the number of timeline samples, /. e. with the dataset size. 
Strictly speaking, the FFT scales log-linearly with time sam- 
ples. However, when performing convolutions, we only take 
the non zero band of N, a tunable but constant factor (see 
Natol i et a/. 200 lb . For the B2K test case under consideration, 
we find that retaining a noise bandwidth of size ~ 10^ is opti- 
mal. We find that about 10^ PCG iterations are needed to reach 
convergence to better than 10"^ precision in the cases under 
consideration. This can be achieved, for instance, in about 5 
minutes on a 128 processor job of an IBM SP3, for the full (8 
PSB, 4 X 10^ total samples) dataset (see also Fig. 

5. Summary and Conclusions 

We have presented our state of the art map-making code to 
jointly reduce multichannel CMB anisotropy and polarisation 
data. ROMA is a parallel (MPI) implementation of the GLS ap- 
proach to map-making, brought to solution by using a PCG it- 
erative method. The only assumptions are that the optical beam 
is purely scalar and axisymmetric, and that the timeline noise 
is (at least piecewise) stationary and uncorrelated across differ- 




Figure 3. Top row: I,Q and U maps obtained by running ROMA on eight "signal-only" simulated timelines (one for each B2K 
bolometer). The display area is restricted to the "shallow" region. Bottom row: difference maps, computed by subtracting the 
above maps from the simulated sky that was input to the simulations. In doing so, the original maps have been degraded to 
NSIDE=1024. Note the change in the colour scale. 
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Figure 4. From left to right, I, Q and U output ROMA maps on eight "signal plus noise" simulated timelines. The assumptions 
for the noise properties are explained in the text. 



ent detectors. For the rest, great care has been taken in tackling 
real world issues, including cross-polarisation, multi-detector 
noise estimation and the problem of missed data. As a test case 
we have reduced with ROMA eight PSB (145 GHz) timelines 
of highly realistic simulated B2K data. We show that the IQU 
maps can be recovered with great precision in the signal-only 
case, while attaining the usual GLS ("optimal") noise suppres- 
sion in the noisy case. We stress that, to our knowledge, this 
is the first joint temperature and polarisation map-making code 
demonstrated to work on a realistic dataset. The code scales lin- 
early with the dataset size, while its parallel behaviour is quasi- 
optimal. It thus represent a viable option to reduce present and 
forthcoming large datasets, including Planck. 
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Figures. Left panel: Difference between the input, simulated sky and the map output by ROMA when processing a single 
detector, "signal-only" (/. e. noiseless) timeline. Both maps (input and output) are intensity maps; the input map was degraded to 
NSIDE=1024 before subtraction. Not surprisingly, the difference is dominated by the polarisation pattern, that is present in the 
timeline but cannot be properly reduced when solving for a single detector map (that is, no polarisation map can be found). Right 
panel: same as left, but here the / map is estimated jointly from all eight detectors {Q and U maps are also computed but are not 
shown here). This figure is de facto the one displayed in left panel, bottom row of Fig. but on a different scale. 
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Figure 6. Shown is inverse time for a single PCG iteration as 
a function of the number of processor elements (A^^^) used by 
ROMA, for the full B2K dataset (eight timelines, consisting of 
5 X 10^ samples each). The machine used is an IBM PowerS 
RS/6000. ROMA scales optimally until N^e ~ 80. 
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